Pylops - Marchenko redatuming

Author: M.Ravasi

In this notebook we consider to extensions of the Marchenko inversion in the case of band-limitied reflectivity (i.e., convolved with wavelet).

We modify the Marchenko equations to account for it and we also use a broad-band (i.e., spiky) initial focusing function.

$$ \begin{bmatrix} \Theta \mathbf{P} \mathbf{f_d^+} \\ \mathbf{0} \end{bmatrix} = \begin{bmatrix} \mathbf{W} & \Theta \mathbf{P} \\ \Theta \mathbf{P^*} & \mathbf{W^*} \end{bmatrix} \begin{bmatrix} \mathbf{f^-} \\ \mathbf{f_m^+} \end{bmatrix} $$

We also consider now the alternating problem of estimating both the focusing functions and the wavelet, using the following equation to update the wavelet.

$$ \begin{bmatrix} \Theta \mathbf{P} \mathbf{f^+} \\ \Theta \mathbf{P^*} \mathbf{f^-} \end{bmatrix} = \begin{bmatrix} \mathbf{F^-} \\ \mathbf{F_m^+} \mathbf{S} \end{bmatrix} \mathbf{w} \rightarrow \mathbf{d_w} = \mathbf{F_w} \mathbf{w} $$

where $\mathbf{P}$ is the bandlimited reflection response (i.e. $P(t, r, s) = w(t) * P(t, r, s)$), $\mathbf{F^-}$ and $\mathbf{F_m^+}$ are convolution operators and $\mathbf{S}$ is a flipping operator that flips the time axis.

In [1]:
%load_ext memory_profiler
%load_ext autoreload
%autoreload 2
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp

from scipy.sparse import csr_matrix, vstack
from scipy.linalg import lstsq, solve
from scipy.sparse.linalg import cg, lsqr
from scipy.signal import convolve, fftconvolve, filtfilt, hilbert

from pylops                            import LinearOperator
from pylops.utils                      import dottest
from pylops.utils.wavelets             import *
from pylops.utils.seismicevents        import *
from pylops.utils.tapers               import *
from pylops.basicoperators             import *
from pylops.signalprocessing           import *
from pylops.waveeqprocessing.mdd       import *
from pylops.waveeqprocessing.marchenko import *
from pylops.optimization.leastsquares  import *
from pylops.optimization.sparsity      import *

Inputs

Input parameters

In [2]:
inputfile = '../data/marchenko/input.npz' # choose file in testdata folder of repo

vel = 2400.0        # velocity
toff = 0.045        # direct arrival time shift
nsmooth = 10        # time window smoothing 
nfmax = 500         # max frequency for MDC (#samples)
nstaper = 11        # source/receiver taper lenght
n_iter = 20         # iterations

jr = 1              # subsampling in r
js = 1              # subsampling in s

clip = 5e5          # figure clipping

Load input

In [3]:
inputdata = np.load(inputfile)

Read and visualize geometry

In [4]:
# Receivers
r = inputdata['r'][:,::jr]
nr = r.shape[1]
dr = r[0,1]-r[0,0]

# Sources
s = inputdata['s'][:,::js]
ns = s.shape[1]
ds = s[0,1]-s[0,0]

# Virtual points
vs = inputdata['vs']

# Density model
rho = inputdata['rho']
z, x = inputdata['z'], inputdata['x']

plt.figure(figsize=(18,9))
plt.imshow(rho, cmap='gray', extent = (x[0], x[-1], z[-1], z[0]))
plt.scatter(s[0, 5::10], s[1, 5::10], marker='*', s=150, c='r', edgecolors='k')
plt.scatter(r[0, ::10],  r[1, ::10], marker='v', s=150, c='b', edgecolors='k')
plt.scatter(vs[0], vs[1], marker='.', s=250, c='m', edgecolors='k')
plt.axis('tight')
plt.xlabel('x [m]'),plt.ylabel('y [m]'),plt.title('Model and Geometry')
plt.xlim(x[0], x[-1]);

Read data

In [5]:
# time axis
t = inputdata['t']
t2 = np.concatenate([-t[::-1], t[1:]])
ot, dt, nt, nt2 = t[0], t[1]-t[0], len(t), len(t2)

# data
R = inputdata['R'][::js, ::jr]
R = np.swapaxes(R, 0, 1) # R[s, r, f] assume particle velocity receiver and integrate over them (via reciprocity)

# tapering
taper = taper3d(nt, [ns, nr], [nstaper, nstaper], tapertype='hanning')
R = R*taper

Define wavelet

In [6]:
# wavelet
ntwav = 21
f0 = 20
phi = 0
amp = .7
ampf = 0.

# original (simple ricker)
wavsymm, twav, wav_c = ricker(t[:ntwav], f0)

# modified wavelet
wav = amp * wavsymm 

WAV = np.fft.rfft(wav, 2**10)
WAVabs = np.abs(WAV)
WAVabsmod = WAVabs + filtfilt(np.ones(10)/10., 1, np.random.normal(0., ampf, len(WAV)), method='pad')
WAV = WAVabsmod * np.exp(1j*np.angle(WAV))
wav = np.fft.irfft(WAV, 2**10)[:len(wav)]

if phi != 0:
    wav = amp * wavsymm(np.cos(np.deg2rad(phi))*wavsymm + np.sin(np.deg2rad(phi))*np.imag(hilbert(wavsymm)))

plt.figure(figsize=(8, 3))
plt.plot(WAVabs, 'k')
plt.plot(WAVabsmod, '--r')
plt.title('Wavelet');
    
plt.figure(figsize=(8, 3))
plt.plot(wavsymm, 'k', label='Original')
plt.plot(wav, 'r', label='Modified')
plt.legend()
plt.title('Wavelet');

Apply wavelet to R

In [7]:
# Convolve R with wavelet
#Rwav = np.apply_along_axis(convolve, 2, R, wav, mode='full')
#Rwav = Rwav[:,:,wav_c:][:, :, :nt]
Rwavsymm = fftconvolve(R, wavsymm[np.newaxis,np.newaxis,:], mode='same', axes=-1)
Rwav = fftconvolve(R, wav[np.newaxis,np.newaxis,:], mode='same', axes=-1)

fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 9))
ax.imshow(R[ns//2].T, cmap='gray', vmin=-1e-2, vmax=1e-2, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax.set_title('R shot=%d' %(ns//2)), ax.set_xlabel(r'$x_R$'), ax.set_ylabel(r'$t$')
ax.axis('tight')
ax.set_ylim(1.5, 0)

fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 9))
ax.imshow(Rwavsymm[ns//2].T, cmap='gray', vmin=-1e-1, vmax=1e-1, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax.set_title('Rwavsymm shot=%d' %(ns//2)), ax.set_xlabel(r'$x_R$'), ax.set_ylabel(r'$t$')
ax.axis('tight')
ax.set_ylim(1.5, 0)

fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 9))
ax.imshow(Rwav[ns//2].T, cmap='gray', vmin=-1e-1, vmax=1e-1, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax.set_title('Rwav shot=%d' %(ns//2)), ax.set_xlabel(r'$x_R$'), ax.set_ylabel(r'$t$')
ax.axis('tight');
ax.set_ylim(1.5, 0);

Read subsurface fields

In [8]:
Gsub = inputdata['Gsub'][:, ::jr]
G0sub = inputdata['G0sub'][:, ::jr]

# convolve G with wavelet
Gsub = np.apply_along_axis(convolve, 0, Gsub, wav, mode='full')
Gsub = Gsub[wav_c:][:nt]
In [9]:
fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub/Gsub.max(), cmap='gray', vmin=-1, vmax=1, 
           extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax2.imshow(G0sub/G0sub.max(), cmap='gray', vmin=-1, vmax=1, 
           extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax2.set_title(r'$G0$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')

ax3.plot(Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(G0sub[:, nr//2]/G0sub.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);
In [10]:
# Direct arrival window - traveltime
toffw = toff - 0.015
directVS = np.sqrt((vs[0]-r[0])**2+(vs[1]-r[1])**2)/vel
directVS_off = directVS - toffw

# Window
idirectVS_off = np.round(directVS_off/dt).astype(np.int)
w = np.zeros((nr, nt))
for ir in range(nr):
    w[ir, :idirectVS_off[ir]]=1            
w = np.hstack((np.fliplr(w), w[:, 1:]))

if nsmooth>0:
    smooth=np.ones(nsmooth)/nsmooth
    w  = filtfilt(smooth, 1, w)    
    
fig, ax = plt.subplots(1, 1,  sharey=True, figsize=(5, 9))
im = ax.imshow(w.T, cmap='gray', extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax.plot(r[0], directVS,'b'),ax.plot(r[0], -directVS,'b')
ax.set_title('Window'), ax.set_xlabel(r'$x_R$'), ax.set_ylabel(r'$t$')
ax.axis('tight')
fig.colorbar(im, ax=ax);

Inversion of Marchenko equations

Use same wavelet for Wop and R

Inversion with modified wavelet

In [11]:
# Add negative time to operators
Rtwosided = np.concatenate((np.zeros((ns, nr, nt-1)), Rwav), axis=-1)
Rtwosided_fft = np.fft.rfft(Rtwosided, 2*nt-1, axis=-1)/np.sqrt(2*nt-1)
Rtwosided_fft = Rtwosided_fft[...,:nfmax]
Rtwosided_fft = Rtwosided_fft.transpose(2, 0, 1) # R[f, s, r]

# Operators
Rop = MDC(Rtwosided_fft, nt=2*nt-1, nv=1, dt=dt, dr=dr, 
          twosided=True, transpose=False, dtype='complex64')
R1op = MDC(Rtwosided_fft, nt=2*nt-1, nv=1, dt=dt, dr=dr, 
           twosided=True, transpose=False, conj=True, dtype='complex64')
Rollop = Roll((2*nt-1) * ns,
              dims=(2*nt-1, ns),
              dir=0, shift=-1)
Wavop = Convolve1D(ns*nt2, wav, wav_c, dims=(nt2, ns), dir=0, dtype='float64')
Wav1op = Convolve1D(ns*nt2, np.flipud(wav), wav_c, dims=(nt2, ns), dir=0, dtype='float64')

Wop = Diagonal(w.T.flatten())

# Input focusing function
fd_plus =  np.concatenate((np.fliplr(G0sub.T).T, np.zeros((nt-1, nr))))
In [12]:
Mop = VStack([HStack([Wavop, -1*Wop*Rop]),
              HStack([-1*Wop*Rollop*R1op, Wav1op])])*BlockDiag([Wop, Wop])
Gop = VStack([HStack([Wavop, -1*Rop]),
              HStack([-1*Rollop*R1op, Wav1op])])

dottest(Gop, 2*ns*(2*nt-1), 2*nr*(2*nt-1), verb=True)
dottest(Mop, 2*ns*(2*nt-1), 2*nr*(2*nt-1), verb=True);
Dot test passed, v^T(Opu)=-1009.002618 - u^T(Op^Tv)=-1009.002618
Dot test passed, v^T(Opu)=280.340654 - u^T(Op^Tv)=280.340654
In [13]:
p0_minus = Rop * fd_plus.flatten()
p0_minus = p0_minus.reshape(nt2, ns).T

d = Wop*Rop*fd_plus.flatten()
d = np.concatenate((d.reshape(nt2, ns), np.zeros((nt2, ns))))

f1_adj = Mop.H*d.flatten()
f1_inv = lsqr(Mop, d.ravel(), iter_lim=n_iter*2, show=True)[0]
res = d.ravel() - Mop * f1_inv

f1_adj = f1_adj.reshape(2*nt2, nr)
f1_inv = f1_inv.reshape(2*nt2, nr)
res = res.reshape(2*nt2, nr)

tau = np.linalg.norm(f1_inv.ravel(), 1)
print('|f1_inv|_1=', tau)

f1_adj_tot = f1_adj + np.concatenate((np.zeros((nt2, nr)), fd_plus))
f1_inv_tot = f1_inv + np.concatenate((np.zeros((nt2, nr)), fd_plus))
 
LSQR            Least-squares solution of  Ax = b
The matrix A has   322998 rows  and   322998 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       40
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.091e+07  2.091e+07    1.0e+00  2.4e-07
     1  0.00000e+00   1.197e+07  1.197e+07    5.7e-01  7.9e-01   6.1e+00  1.0e+00
     2  0.00000e+00   8.305e+06  8.305e+06    4.0e-01  3.8e-01   9.7e+00  2.3e+00
     3  0.00000e+00   6.502e+06  6.502e+06    3.1e-01  2.9e-01   1.2e+01  3.7e+00
     4  0.00000e+00   5.145e+06  5.145e+06    2.5e-01  2.2e-01   1.4e+01  5.4e+00
     5  0.00000e+00   4.316e+06  4.316e+06    2.1e-01  1.8e-01   1.6e+01  7.1e+00
     6  0.00000e+00   3.620e+06  3.620e+06    1.7e-01  1.8e-01   1.7e+01  9.0e+00
     7  0.00000e+00   3.097e+06  3.097e+06    1.5e-01  1.3e-01   1.8e+01  1.1e+01
     8  0.00000e+00   2.783e+06  2.783e+06    1.3e-01  1.1e-01   2.0e+01  1.3e+01
     9  0.00000e+00   2.511e+06  2.511e+06    1.2e-01  1.1e-01   2.1e+01  1.5e+01
    10  0.00000e+00   2.273e+06  2.273e+06    1.1e-01  9.7e-02   2.2e+01  1.8e+01
    30  0.00000e+00   9.389e+05  9.389e+05    4.5e-02  2.8e-02   3.9e+01  8.4e+01
    31  0.00000e+00   9.189e+05  9.189e+05    4.4e-02  2.6e-02   3.9e+01  8.9e+01
    32  0.00000e+00   8.994e+05  8.994e+05    4.3e-02  2.5e-02   4.0e+01  9.3e+01
    33  0.00000e+00   8.821e+05  8.821e+05    4.2e-02  2.5e-02   4.1e+01  9.7e+01
    34  0.00000e+00   8.646e+05  8.646e+05    4.1e-02  2.4e-02   4.1e+01  1.0e+02
    35  0.00000e+00   8.492e+05  8.492e+05    4.1e-02  2.2e-02   4.2e+01  1.1e+02
    36  0.00000e+00   8.354e+05  8.354e+05    4.0e-02  2.2e-02   4.2e+01  1.1e+02
    37  0.00000e+00   8.219e+05  8.219e+05    3.9e-02  2.1e-02   4.3e+01  1.2e+02
    38  0.00000e+00   8.094e+05  8.094e+05    3.9e-02  2.0e-02   4.4e+01  1.2e+02
    39  0.00000e+00   7.984e+05  7.984e+05    3.8e-02  1.8e-02   4.4e+01  1.2e+02
    40  0.00000e+00   7.874e+05  7.874e+05    3.8e-02  1.8e-02   4.5e+01  1.3e+02
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 7.9e+05   anorm = 4.5e+01   arnorm = 6.5e+05
itn   =      40   r2norm = 7.9e+05   acond = 1.3e+02   xnorm  = 7.0e+06
 
|f1_inv|_1= 664815486.2219818
In [14]:
g_adj = Gop*f1_adj_tot.flatten()
g_inv = Gop*f1_inv_tot.flatten()

g_adj = g_adj.reshape(2*nt2, ns)
g_inv = g_inv.reshape(2*nt2, ns)

f1_adj_minus, f1_adj_plus =  f1_adj_tot[:nt2].T, f1_adj_tot[nt2:].T
f1_inv_minus, f1_inv_plus =  f1_inv_tot[:nt2].T, f1_inv_tot[nt2:].T

g_adj_minus, g_adj_plus =  -g_adj[:nt2].T, np.fliplr(g_adj[nt2:].T)
g_inv_minus, g_inv_plus =  -g_inv[:nt2].T, np.fliplr(g_inv[nt2:].T)

g_inv_tot = g_inv_minus + g_inv_plus
In [15]:
# Need to recreate combined data as new implementation stacks over time instead of space
d_plot = np.concatenate((d[:(2*nt-1)], d[(2*nt-1):]), axis=1).T
res_plot = np.concatenate((res[:(2*nt-1)], res[(2*nt-1):]), axis=1).T
f1_adj_tot_plot = np.concatenate((f1_adj_tot[:(2*nt-1)], f1_adj_tot[(2*nt-1):]), axis=1).T
f1_inv_tot_plot = np.concatenate((f1_inv_tot[:(2*nt-1)], f1_inv_tot[(2*nt-1):]), axis=1).T
g_adj_plot = np.concatenate((g_adj[:(2*nt-1)], g_adj[(2*nt-1):]), axis=1).T
g_inv_plot = np.concatenate((g_inv[:(2*nt-1)], g_inv[(2*nt-1):]), axis=1).T

fig, axs = plt.subplots(1, 6, sharey=True, figsize=(16, 7))
axs[0].imshow(d_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[0].set_title('d'), axs[0].set_xlabel(r'$x_R$'), axs[3].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1)
axs[1].imshow(res_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[1].set_title('res'), axs[0].set_xlabel(r'$x_R$'), axs[3].set_ylabel(r'$t$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1)
axs[2].imshow(f1_adj_tot_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f_{adj}$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
axs[3].imshow(f1_inv_tot_plot.T, cmap='gray', vmin=-clip/10, vmax=clip/10, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[3].set_title(r'$f_{inv}$'), axs[0].set_xlabel(r'$x_R$')
axs[3].axis('tight')
axs[3].set_ylim(1, -1);
axs[4].imshow(g_adj_plot.T, cmap='gray', vmin=-10*clip, vmax=10*clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[4].set_title(r'$g_{adj}$'), axs[0].set_xlabel(r'$x_R$')
axs[4].axis('tight')
axs[4].set_ylim(1, -1);
axs[5].imshow(g_inv_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[5].set_title(r'$g_{inv}$'), axs[0].set_xlabel(r'$x_R$')
axs[5].axis('tight')
axs[5].set_ylim(1, -1);
In [16]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_adj_minus.T, cmap='gray', vmin=-10*clip, vmax=10*clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$ adj'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_adj_plus.T, cmap='gray', vmin=-10*clip, vmax=10*clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$ adj'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_inv_minus.T, cmap='gray', vmin=-clip/10, vmax=clip/10, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_inv_plus.T, cmap='gray', vmin=-clip/10, vmax=clip/10, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
In [17]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow((1-w.T)*g_inv_minus.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_inv_plus.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);

fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-clip, vmax=clip, 
           extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot.T, cmap='gray', vmin=-clip, vmax=clip,
           extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(t**2*Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(t**2*g_inv_tot[nr//2, nt-1:]/g_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

Use different wavelet for Wop and R

Let's what happens if we use the wrong wavelet (use wavsymm for the wavelet operator but R with the modified wavelet)

In [18]:
Wavsymmop = Convolve1D(ns*nt2, wavsymm, wav_c, dims=(nt2, ns), 
                       dir=0, dtype='float64')
Wavsymm1op = Convolve1D(ns*nt2, np.flipud(wavsymm), wav_c, dims=(nt2, ns), 
                        dir=0, dtype='float64')

# Input focusing function
fd_plus =  np.concatenate((np.fliplr(G0sub.T).T, np.zeros((nt-1, nr))))

Msymmop = VStack([HStack([Wavsymmop, -1*Wop*Rop]),
                  HStack([-1*Wop*Rollop*R1op, Wavsymm1op])])*BlockDiag([Wop, Wop])
Gsymmop = VStack([HStack([Wavsymmop, -1*Rop]),
                  HStack([-1*Rollop*R1op, Wavsymm1op])])

p0_symm_minus = Rop * fd_plus.flatten()
p0_symm_minus = p0_symm_minus.reshape(nt2, ns).T

d_symm = Wop * Rop * fd_plus.flatten()
d_symm = np.concatenate((d_symm.reshape(nt2, ns), np.zeros((nt2, ns))))

f1_symm_adj = Msymmop.H*d_symm.flatten()
f1_symm_inv = lsqr(Msymmop, d_symm.ravel(), iter_lim=n_iter*2, show=True)[0]
res_symm = d.ravel() - Msymmop * f1_symm_inv

f1_symm_adj = f1_symm_adj.reshape(2*nt2, nr)
f1_symm_inv = f1_symm_inv.reshape(2*nt2, nr)
res_symm = res_symm.reshape(2*nt2, nr)

f1_symm_adj_tot = f1_symm_adj + np.concatenate((np.zeros((nt2, nr)), fd_plus))
f1_symm_inv_tot = f1_symm_inv + np.concatenate((np.zeros((nt2, nr)), fd_plus))

g_symm_adj = Gsymmop*f1_symm_adj_tot.flatten()
g_symm_inv = Gsymmop*f1_symm_inv_tot.flatten()

g_symm_adj = g_symm_adj.reshape(2*nt2, ns)
g_symm_inv = g_symm_inv.reshape(2*nt2, ns)

f1_symm_adj_minus, f1_symm_adj_plus =  f1_symm_adj_tot[:nt2].T, f1_symm_adj_tot[nt2:].T
f1_symm_inv_minus, f1_symm_inv_plus =  f1_symm_inv_tot[:nt2].T, f1_symm_inv_tot[nt2:].T

g_symm_adj_minus, g_symm_adj_plus =  -g_symm_adj[:nt2].T, np.fliplr(g_symm_adj[nt2:].T)
g_symm_inv_minus, g_symm_inv_plus =  -g_symm_inv[:nt2].T, np.fliplr(g_symm_inv[nt2:].T)

g_symm_inv_tot = g_symm_inv_minus + g_symm_inv_plus
 
LSQR            Least-squares solution of  Ax = b
The matrix A has   322998 rows  and   322998 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       40
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.091e+07  2.091e+07    1.0e+00  3.4e-07
     1  0.00000e+00   1.045e+07  1.045e+07    5.0e-01  7.9e-01   8.1e+00  1.0e+00
     2  0.00000e+00   6.975e+06  6.975e+06    3.3e-01  4.1e-01   1.2e+01  2.2e+00
     3  0.00000e+00   5.077e+06  5.077e+06    2.4e-01  3.1e-01   1.5e+01  3.6e+00
     4  0.00000e+00   4.029e+06  4.029e+06    1.9e-01  2.2e-01   1.8e+01  5.1e+00
     5  0.00000e+00   3.315e+06  3.315e+06    1.6e-01  1.9e-01   2.0e+01  6.8e+00
     6  0.00000e+00   2.768e+06  2.768e+06    1.3e-01  1.5e-01   2.1e+01  8.7e+00
     7  0.00000e+00   2.450e+06  2.450e+06    1.2e-01  1.4e-01   2.3e+01  1.1e+01
     8  0.00000e+00   2.149e+06  2.149e+06    1.0e-01  1.2e-01   2.5e+01  1.3e+01
     9  0.00000e+00   1.936e+06  1.936e+06    9.3e-02  1.0e-01   2.6e+01  1.5e+01
    10  0.00000e+00   1.756e+06  1.756e+06    8.4e-02  9.8e-02   2.8e+01  1.7e+01
    30  0.00000e+00   7.040e+05  7.040e+05    3.4e-02  2.8e-02   4.8e+01  8.4e+01
    31  0.00000e+00   6.873e+05  6.873e+05    3.3e-02  2.7e-02   4.9e+01  8.9e+01
    32  0.00000e+00   6.726e+05  6.726e+05    3.2e-02  2.6e-02   5.0e+01  9.3e+01
    33  0.00000e+00   6.587e+05  6.587e+05    3.1e-02  2.5e-02   5.1e+01  9.7e+01
    34  0.00000e+00   6.457e+05  6.457e+05    3.1e-02  2.4e-02   5.1e+01  1.0e+02
    35  0.00000e+00   6.330e+05  6.330e+05    3.0e-02  2.3e-02   5.2e+01  1.1e+02
    36  0.00000e+00   6.216e+05  6.216e+05    3.0e-02  2.2e-02   5.3e+01  1.1e+02
    37  0.00000e+00   6.109e+05  6.109e+05    2.9e-02  2.1e-02   5.4e+01  1.2e+02
    38  0.00000e+00   6.011e+05  6.011e+05    2.9e-02  2.1e-02   5.4e+01  1.2e+02
    39  0.00000e+00   5.917e+05  5.917e+05    2.8e-02  2.0e-02   5.5e+01  1.2e+02
    40  0.00000e+00   5.825e+05  5.825e+05    2.8e-02  2.0e-02   5.6e+01  1.3e+02
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 5.8e+05   anorm = 5.6e+01   arnorm = 6.5e+05
itn   =      40   r2norm = 5.8e+05   acond = 1.3e+02   xnorm  = 4.6e+06
 
In [19]:
# Need to recreate combined data as new implementation stacks over time instead of space
res_plot = np.concatenate((res_symm[:(2*nt-1)], res_symm[(2*nt-1):]), axis=1).T
f1_adj_tot_plot = np.concatenate((f1_symm_adj_tot[:(2*nt-1)], f1_symm_adj_tot[(2*nt-1):]), axis=1).T
f1_inv_tot_plot = np.concatenate((f1_symm_inv_tot[:(2*nt-1)], f1_symm_inv_tot[(2*nt-1):]), axis=1).T
g_adj_plot = np.concatenate((g_symm_adj[:(2*nt-1)], g_symm_adj[(2*nt-1):]), axis=1).T
g_inv_plot = np.concatenate((g_symm_inv[:(2*nt-1)], g_symm_inv[(2*nt-1):]), axis=1).T

fig, axs = plt.subplots(1, 6, sharey=True, figsize=(16, 7))
axs[0].imshow(d_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[0].set_title('d'), axs[0].set_xlabel(r'$x_R$'), axs[3].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1)
axs[1].imshow(res_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[1].set_title('res'), axs[0].set_xlabel(r'$x_R$'), axs[3].set_ylabel(r'$t$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1)
axs[2].imshow(f1_adj_tot_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f_{adj}$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
axs[3].imshow(f1_inv_tot_plot.T, cmap='gray', vmin=-clip/10, vmax=clip/10, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[3].set_title(r'$f_{inv}$'), axs[0].set_xlabel(r'$x_R$')
axs[3].axis('tight')
axs[3].set_ylim(1, -1);
axs[4].imshow(g_adj_plot.T, cmap='gray', vmin=-10*clip, vmax=10*clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[4].set_title(r'$g_{adj}$'), axs[0].set_xlabel(r'$x_R$')
axs[4].axis('tight')
axs[4].set_ylim(1, -1);
axs[5].imshow(g_inv_plot.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[5].set_title(r'$g_{inv}$'), axs[0].set_xlabel(r'$x_R$')
axs[5].axis('tight')
axs[5].set_ylim(1, -1);
In [20]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_symm_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_symm_adj_minus.T, cmap='gray', vmin=-10*clip, vmax=10*clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$ adj'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_symm_adj_plus.T, cmap='gray', vmin=-10*clip, vmax=10*clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$ adj'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_symm_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_symm_inv_minus.T, cmap='gray', vmin=-clip/10, vmax=clip/10, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_symm_inv_plus.T, cmap='gray', vmin=-clip/10, vmax=clip/10, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
In [21]:
fig, ax = plt.subplots(1, 1, figsize=(12, 2))
ax.plot(Rwavsymm[50, 50], 'k', label='Original')
ax.plot(Rwav[50, 50], 'r', label='Modified')
ax.set_title('R')
ax.set_xlim(0, 300)
ax.legend()

fig, axs = plt.subplots(2, 1, figsize=(12, 5))
axs[0].plot(f1_symm_inv_minus[50], 'k')
axs[0].plot(f1_inv_minus[50], 'r')
axs[0].set_title('f^-')
axs[0].set_xlim(600, 900)
axs[1].plot(f1_symm_inv_plus[50], 'k', label='Original')
axs[1].plot(f1_inv_plus[50], 'r', label='Modified')
axs[1].set_title('f^+')
axs[1].set_xlim(600, 900)
axs[1].legend();
In [22]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_symm_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow((1-w.T)*g_symm_inv_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_symm_inv_plus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);

fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_symm_inv_tot.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(t**2*Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(t**2*g_symm_inv_tot[nr//2, nt-1:]/g_symm_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

We compare now the results

In [23]:
fig, axs = plt.subplots(1, 2, sharey=True, figsize=(10, 9))
axs[0].imshow(f1_adj_minus.T - f1_symm_inv_minus.T, cmap='gray', 
              vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$f^-$ error adj'), axs[0].set_xlabel(r'$x_R$')
axs[0].plot(r[0], directVS,'b'),axs[0].plot(r[0], -directVS,'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_adj_plus.T - f1_symm_adj_plus.T, cmap='gray', 
              vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^+$ error adj'), axs[1].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
In [24]:
fig, axs = plt.subplots(1, 2, sharey=True, figsize=(10, 9))
axs[0].imshow(f1_inv_minus.T - f1_symm_inv_minus.T, cmap='gray', 
              vmin=-clip/10, vmax=clip/10, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$f^-$ error inv'), axs[0].set_xlabel(r'$x_R$')
axs[0].plot(r[0], directVS,'b'),axs[0].plot(r[0], -directVS,'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_inv_plus.T - f1_symm_inv_plus.T, cmap='gray', 
              vmin=-clip/10, vmax=clip/10, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^+$ error inv'), axs[1].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
In [25]:
fig, axs = plt.subplots(1, 2, sharey=True, figsize=(10, 9))
axs[0].imshow(g_inv_minus.T - g_symm_inv_minus.T, cmap='gray', 
              vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$g^-$ error inv'), axs[0].set_xlabel(r'$x_R$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow(g_inv_plus.T + g_symm_inv_plus.T, cmap='gray', 
              vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^+$ error inv'), axs[1].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);

We can see how in both cases we have been able to come up with a set of focusing functions that can explain the data (near to zero residual).

However only the first solution is correct and leads to no artefacts in $G^-$ while the second does. Althought not always valid, we can perhaps try to use this argument (minimum L2 norm of $G^-$) as a regularization term when inverting for the wavelet?

Wavelet estimation equations

Let's try to invert for the wavelet given the already estimated focusing functions

In [26]:
f1_inv_tot = f1_inv + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus))
f1_inv_minus, f1_inv_plus =  f1_inv_tot[:nt2].T, f1_inv_tot[nt2:].T

d1 = Wop*Rop*f1_inv_plus.T.flatten()
d2 = Wop*R1op*f1_inv_minus.T.flatten()
dw = np.concatenate((d1.reshape(nt2, ns).T, d2.reshape(nt2, ns).T))
In [27]:
W1op = Diagonal(w.flatten())
F1plus = VStack([Convolve1D(nt2, f1_inv[nt2:, ir], offset=nt-2) for ir in range(nr)])
F1minus = VStack([Convolve1D(nt2, f1_inv[:nt2, ir], offset=nt-1) for ir in range(nr)])
S = Flip(nt2)
Fw = VStack([W1op*F1minus, W1op*F1plus*S])
In [28]:
extended_wav = np.pad(wav, (nt-ntwav, nt-ntwav), mode='constant')

dwpred = Fw * extended_wav
dwpred = dwpred.reshape(2 * nr, nt2)
In [29]:
fig, axs = plt.subplots(1, 3, figsize=(15, 9))
axs[0].imshow(dw.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[0].axis('tight');
axs[1].imshow(dwpred.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[1].axis('tight')
axs[2].imshow(dw.T - dwpred.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[2].axis('tight');

plt.figure(figsize=(12, 3))
plt.plot(dw[50])
plt.plot(dwpred[50]);
plt.xlim(600, 800);
plt.figure(figsize=(12, 3))
plt.plot(dw[150])
plt.plot(dwpred[150]);
plt.xlim(600, 800);
In [30]:
wavest = lsqr(Fw, dw.flatten(), iter_lim=5, show=True)[0]
 
LSQR            Least-squares solution of  Ax = b
The matrix A has   322998 rows  and     1599 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =        5
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.509e+07  2.509e+07    1.0e+00  6.1e-01
     1  0.00000e+00   6.721e+06  6.721e+06    2.7e-01  8.3e-01   1.6e+07  1.0e+00
     2  0.00000e+00   2.876e+06  2.876e+06    1.1e-01  5.0e-01   2.2e+07  2.1e+00
     3  0.00000e+00   1.661e+06  1.661e+06    6.6e-02  2.7e-01   2.6e+07  3.2e+00
     4  0.00000e+00   1.395e+06  1.395e+06    5.6e-02  1.9e-01   2.9e+07  4.4e+00
     5  0.00000e+00   1.227e+06  1.227e+06    4.9e-02  1.0e-01   3.3e+07  6.0e+00
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 1.2e+06   anorm = 3.3e+07   arnorm = 4.0e+12
itn   =       5   r2norm = 1.2e+06   acond = 6.0e+00   xnorm  = 1.7e+00
 
In [31]:
plt.figure()
plt.plot(twav, wavsymm, '--k', lw=1, label='Modified (used)')
plt.plot(twav, extended_wav[nt-ntwav:nt+ntwav-1], 'k', lw=2, label='Original')
plt.plot(twav, wavest[nt-ntwav:nt+ntwav-1], '--r', label='Estimated')
plt.grid()

You can see that the correct wavelet has been estimated.

We now use the focusing functions retrieved using the the wrong wavelet (which could have been the initial guess in the alternating algorithm)

In [32]:
f1_symm_inv_tot = f1_symm_inv + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus))
f1_symm_inv_minus, f1_symm_inv_plus =  f1_symm_inv_tot[:nt2].T, f1_symm_inv_tot[nt2:].T

d1_symm = Wop*Rop*f1_symm_inv_plus.T.flatten()
d2_symm = Wop*R1op*f1_symm_inv_minus.T.flatten()
dw_symm = np.concatenate((d1_symm.reshape(nt2, ns).T, d2_symm.reshape(nt2, ns).T))

W1op = Diagonal(w.flatten())
F1plus = VStack([Convolve1D(nt2, f1_symm_inv[nt2:, ir], offset=nt-2) for ir in range(nr)])
F1minus = VStack([Convolve1D(nt2, f1_symm_inv[:nt2, ir], offset=nt-1) for ir in range(nr)])
S = Flip(nt2)
Fw = VStack([W1op*F1minus, W1op*F1plus*S])

extended_wav = np.pad(wavsymm, (nt-ntwav, nt-ntwav), mode='constant')

dwpred_symm = Fw * extended_wav
dwpred_symm = dwpred_symm.reshape(2 * nr, nt2)
In [33]:
fig, axs = plt.subplots(1, 3, figsize=(15, 9))
axs[0].imshow(dw_symm.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[0].axis('tight');
axs[1].imshow(dwpred_symm.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[1].axis('tight')
axs[2].imshow(dw_symm.T - dwpred_symm.T, cmap='gray', vmin=-np.abs(dw).max(), vmax=np.abs(dw).max())
axs[2].axis('tight');

plt.figure(figsize=(12, 3))
plt.plot(dw_symm[50])
plt.plot(dwpred_symm[50]);
plt.xlim(600, 800);
plt.figure(figsize=(12, 3))
plt.plot(dw_symm[150])
plt.plot(dwpred_symm[150]);
plt.xlim(600, 800);
In [34]:
wavest_symm = lsqr(Fw, dw.flatten(), iter_lim=5, show=True)[0]

plt.figure()
plt.plot(twav, wavsymm, '--k', lw=1, label='Modified (used)')
plt.plot(twav, extended_wav[nt-ntwav:nt+ntwav-1], 'k', lw=2, label='Original')
plt.plot(twav, wavest_symm[nt-ntwav:nt+ntwav-1], '--r', label='Estimaged')
plt.grid()
 
LSQR            Least-squares solution of  Ax = b
The matrix A has   322998 rows  and     1599 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =        5
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.509e+07  2.509e+07    1.0e+00  3.8e-01
     1  0.00000e+00   7.629e+06  7.629e+06    3.0e-01  7.5e-01   1.0e+07  1.0e+00
     2  0.00000e+00   4.696e+06  4.696e+06    1.9e-01  3.0e-01   1.4e+07  2.1e+00
     3  0.00000e+00   4.127e+06  4.127e+06    1.6e-01  1.2e-01   1.7e+07  3.2e+00
     4  0.00000e+00   4.014e+06  4.014e+06    1.6e-01  5.9e-02   1.9e+07  4.4e+00
     5  0.00000e+00   3.958e+06  3.958e+06    1.6e-01  3.6e-02   2.1e+07  6.0e+00
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 4.0e+06   anorm = 2.1e+07   arnorm = 3.0e+12
itn   =       5   r2norm = 4.0e+06   acond = 6.0e+00   xnorm  = 2.7e+00
 

We can see that the inversion returns a completely new wavelet when using the wrong focusing functions. This means that if we let the focusing functions converge (zero residual) and estimate the wavelet at that point we will get out something that is compensating for the fact that we have used incosistent wavelets in the data and Wop.

We need to try using a sparse solver for the focusing functions to have slow convergence and use an alternating optimization as in EPSI.

Here we just try the sparse solver and let it converge

In [52]:
# Initial focusing function
iG0 = np.argmax(G0sub, axis=0)
G0sub_sparse = np.zeros_like(G0sub)
G0sub_sparse[iG0, np.arange(nr)] = G0sub.max()
G0sub_sparse = fftconvolve(G0sub_sparse, np.hanning(7)[:, np.newaxis], mode='same', axes=0)
#G0sub_sparse = G0sub.copy()

fd_plus_sparse =  np.concatenate((np.fliplr(G0sub_sparse.T).T, np.zeros((nt-1, nr))))
d_sparse = Wop*Rop*fd_plus_sparse.flatten()
d_sparse = np.concatenate((d_sparse.reshape(nt2, ns), np.zeros((nt2, ns))))
In [53]:
f1_sparse_inv, _, info = SPGL1(Mop, d_sparse.flatten(), tau=tau, iter_lim=100, verbosity=2)
f1_sparse_inv = f1_sparse_inv.reshape(2*nt2, nr)
print('|f1_inv|_1=', np.linalg.norm(f1_sparse_inv.ravel(), 1))

f1_sparse_inv_tot = f1_sparse_inv + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus_sparse))
================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 6.65e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 6.65e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :      100

iterr      Objective   Relative Gap      gnorm   stepg   nnz_x   nnz_g

    0  2.8148534e+07  8.9552701e+00   5.34e+06     0.0       0       0       
    1  2.8148264e+07  8.9553391e+00   5.34e+06     0.0   66223       1       
    2  1.5686899e+07  1.1883190e+01   2.20e+06     0.0   86044       0       
    3  1.2923957e+07  1.5801482e+01   2.10e+06     0.0   86094       1       
    4  9.8180576e+06  1.6689291e+01   1.27e+06     0.0   86120       0       
    5  8.3359394e+06  1.1279926e+01   6.93e+05     0.0   86130       1       
    6  7.6632130e+06  1.0389754e+01   5.47e+05     0.0   86135       1       
    7  6.1448168e+06  1.2820593e+01   4.20e+05     0.0   86138       1       
    8  6.2380750e+06  3.9827843e+01   1.18e+06    -0.3   86164       1       
    9  6.2602221e+06  4.0449400e+01   1.22e+06     0.0   86167       0       
   10  4.7457196e+06  1.1396488e+01   2.30e+05     0.0   86173       1       
   11  4.6206815e+06  1.0579222e+01   2.06e+05     0.0   86179       1       
   12  3.8949334e+06  1.1144352e+01   1.52e+05     0.0   86197       1       
   13  3.2092737e+06  6.3609519e+01   5.08e+05     0.0   44202       1       
   14  3.4902472e+06  7.2911188e+01   6.71e+05    -0.9   50562       1       
   15  2.6607652e+06  1.3345511e+01   8.53e+04     0.0   69634       1       
   16  2.6310320e+06  1.0437795e+01   6.77e+04     0.0   60923       1       
   17  2.5727186e+06  9.7426998e+00   6.07e+04     0.0   54262       1       
   18  2.2102280e+06  3.7444628e+01   1.48e+05     0.0   42676       1       
   19  2.4976719e+06  1.0874330e+02   5.21e+05    -0.9   51314       1       
   20  2.1704998e+06  4.9799044e+01   1.87e+05     0.0   56663       1       
   21  2.0897301e+06  1.0994439e+01   4.79e+04     0.0   60001       1       
   22  2.0795121e+06  9.7483442e+00   4.32e+04     0.0   55181       1       
   23  2.0281558e+06  9.2522752e+00   3.92e+04     0.0   45975       1       
   24  1.9490909e+06  6.4457863e+01   1.93e+05    -0.3   39815       1       
   25  2.0498603e+06  8.8242441e+01   2.88e+05    -0.9   48785       0       
   26  1.8360064e+06  1.9548363e+01   5.82e+04     0.0   53507       1       
   27  1.8236903e+06  9.2932659e+00   3.21e+04     0.0   50367       1       
   28  1.8170457e+06  8.8039029e+00   3.06e+04     0.0   46932       1       
   29  1.6844342e+06  2.1731499e+01   5.39e+04     0.0   37566       1       
   30  1.7511299e+06  7.1548793e+01   1.71e+05    -1.2   42896       1       
   31  1.6903682e+06  4.7425197e+01   1.10e+05     0.0   50433       1       
   32  1.6493034e+06  8.7864976e+00   2.56e+04     0.0   50948       1       
   33  1.6459424e+06  8.5177275e+00   2.49e+04     0.0   47641       1       
   34  1.6276583e+06  8.5986383e+00   2.43e+04     0.0   41108       1       
   35  1.6295643e+06  1.0514773e+02   2.18e+05     0.0   35226       1       
   36  1.5596828e+06  8.7549410e+01   1.64e+05    -1.2   48036       1       
   37  1.4800168e+06  3.2558318e+01   5.99e+04     0.0   55796       1       
   38  1.4668613e+06  8.7881604e+00   2.00e+04     0.0   52362       1       
   39  1.4642297e+06  7.9678755e+00   1.86e+04     0.0   49276       1       
   40  1.4514630e+06  7.6740553e+00   1.74e+04     0.0   41658       1       
   41  1.4508109e+06  3.2365005e+01   5.73e+04    -0.3   40440       0       
   42  1.4380410e+06  2.1614346e+01   3.89e+04    -0.6   45592       1       
   43  1.4332062e+06  8.5866046e+00   1.87e+04     0.0   44737       1       
   44  1.4314553e+06  7.3203298e+00   1.66e+04     0.0   42377       1       
   45  1.4241971e+06  8.0558999e+00   1.75e+04     0.0   39348       1       
   46  1.4294767e+06  6.1361664e+01   9.88e+04    -0.3   38038       1       
   47  1.4248414e+06  5.3303705e+01   8.70e+04    -0.6   43065       1       
   48  1.4016605e+06  8.0206435e+00   1.70e+04     0.0   45681       1       
   49  1.4004056e+06  7.0496647e+00   1.55e+04     0.0   42444       1       
   50  1.3959643e+06  6.4212090e+00   1.44e+04     0.0   39109       1       
   51  1.3581506e+06  4.3995908e+01   6.60e+04     0.0   35996       1       
   52  1.3616211e+06  6.6218766e+01   9.60e+04    -1.5   41918       1       
   53  1.3460072e+06  1.9042934e+01   3.06e+04     0.0   45537       1       
   54  1.3437170e+06  6.1481305e+00   1.28e+04     0.0   43686       1       
   55  1.3427028e+06  6.1467404e+00   1.27e+04     0.0   41429       1       
   56  1.3336602e+06  9.8453941e+00   1.74e+04     0.0   37507       1       
   57  1.3354747e+06  3.4469181e+01   5.09e+04    -1.2   38235       1       
   58  1.3339349e+06  3.3029251e+01   4.82e+04     0.0   39413       1       
   59  1.3281978e+06  6.6361910e+00   1.31e+04     0.0   41023       1       
   60  1.3274935e+06  5.9761137e+00   1.22e+04     0.0   39381       1       
   61  1.3226107e+06  6.3777164e+00   1.26e+04     0.0   37472       1       
   62  1.3114106e+06  1.0456016e+02   1.38e+05    -0.6   35609       1       
   63  1.2876109e+06  4.1586662e+01   5.59e+04    -1.5   42860       1       
   64  1.2801147e+06  9.4470057e+00   1.53e+04     0.0   47218       1       
   65  1.2792063e+06  5.4344794e+00   1.04e+04     0.0   44011       1       
   66  1.2778125e+06  5.1221083e+00   9.96e+03     0.0   41043       1       
   67  1.2696459e+06  2.0343783e+01   2.85e+04     0.0   37177       1       
   68  1.2707743e+06  2.6617762e+01   3.56e+04    -1.2   39419       0       
   69  1.2665145e+06  9.0523713e+00   1.46e+04     0.0   41136       1       
   70  1.2656993e+06  5.4890521e+00   1.02e+04     0.0   40155       1       
   71  1.2650859e+06  5.6644469e+00   1.04e+04     0.0   39159       1       
   72  1.2453201e+06  2.1613202e+01   2.83e+04     0.0   36213       0       
   73  1.2455482e+06  2.5595496e+01   3.34e+04    -1.8   39580       1       
   74  1.2413527e+06  7.2592601e+00   1.16e+04     0.0   41636       1       
   75  1.2407202e+06  4.8268142e+00   8.87e+03     0.0   40327       1       
   76  1.2400206e+06  4.6270767e+00   8.63e+03     0.0   39101       1       
   77  1.2310094e+06  1.9707282e+01   2.58e+04     0.0   36476       1       
   78  1.2329500e+06  3.8351785e+01   4.68e+04    -1.2   40953       0       
   79  1.2282837e+06  1.8909404e+01   2.48e+04     0.0   43041       1       
   80  1.2261926e+06  5.4900177e+00   9.39e+03     0.0   42086       1       
   81  1.2258375e+06  5.7051872e+00   9.62e+03     0.0   41130       1       
   82  1.2202754e+06  6.3536257e+00   1.01e+04     0.0   37548       1       
   83  1.2203032e+06  2.5543276e+01   3.19e+04    -0.9   38905       1       
   84  1.2187782e+06  1.9179149e+01   2.43e+04    -0.3   40733       1       
   85  1.2172457e+06  4.8873300e+00   8.49e+03     0.0   41293       1       
   86  1.2169441e+06  4.5479887e+00   8.09e+03     0.0   39767       1       
   87  1.2151749e+06  4.9349448e+00   8.50e+03     0.0   37696       1       
   88  1.2107162e+06  4.7164689e+01   5.44e+04    -0.6   36156       1       
   89  1.2168200e+06  5.1379428e+01   6.06e+04    -1.2   40692       1       
   90  1.2028309e+06  6.4477665e+00   9.91e+03     0.0   47590       1       
   91  1.2023660e+06  4.5942548e+00   7.89e+03     0.0   43824       1       
   92  1.2019759e+06  4.5275736e+00   7.78e+03     0.0   40902       1       
   93  1.1945340e+06  1.8977860e+01   2.32e+04     0.0   36081       0       
   94  1.1958311e+06  2.5046383e+01   2.95e+04    -1.5   39226       1       
   95  1.1925141e+06  8.7570746e+00   1.22e+04     0.0   42101       1       
   96  1.1919019e+06  3.7705464e+00   6.76e+03     0.0   40824       1       
   97  1.1916613e+06  3.9118512e+00   6.91e+03     0.0   39710       1       
   98  1.1853029e+06  1.0099936e+01   1.33e+04     0.0   36385       0       
   99  1.1856392e+06  2.1495798e+01   2.55e+04    -1.5   39143       1       
Restoring best iterate to objective 1185302.9020167666

ERROR EXIT -- Too many iterations

Products with A     :      77        Total time   (secs) :    96.1
Products with A^H   :     101        Project time (secs) :     9.9
Newton iterations   :       0        Mat-vec time (secs) :    81.2
Line search its     :      75        Subspace iterations :       0
|f1_inv|_1= 664815486.2219867
In [54]:
plt.figure()
plt.plot(info['xnorm1'])
plt.axhline(tau, color='r');
In [55]:
g_sparse_inv = Gop*f1_sparse_inv_tot.flatten()
g_sparse_inv = g_sparse_inv.reshape(2*nt2, ns)

f1_sparse_inv_minus, f1_sparse_inv_plus =  f1_sparse_inv_tot[:nt2].T, f1_sparse_inv_tot[nt2:].T
g_sparse_inv_minus, g_sparse_inv_plus =  -g_sparse_inv[:nt2].T, np.fliplr(g_sparse_inv[nt2:].T)
g_sparse_inv_tot = g_sparse_inv_minus + g_sparse_inv_plus
In [56]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_sparse_inv_minus.T, cmap='gray', vmin=-clip/10, vmax=clip/10, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_sparse_inv_plus.T, cmap='gray', vmin=-clip/10, vmax=clip/10, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
In [57]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow((1-w.T)*g_sparse_inv_minus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_sparse_inv_plus.T, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);

fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-clip, vmax=clip, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_sparse_inv_tot.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(t**2*Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(t**2*g_sparse_inv_tot[nr//2, nt-1:]/g_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

Alternating optimization

Finally we perform alternating optimization to estimate both focusing functions and wavelet

In [64]:
optimalsigma = 0.01 * np.linalg.norm(d.ravel())
nouter = 10
quickplot=True

# initial wavelet
ntwav = len(twav)//2 + 1
nfft = 2**10
In [65]:
# Estimate wavelet spectrum
wavest0_fft = np.mean(np.abs(np.fft.fft(Rwav, nfft, axis=-1)), axis=(0, 1))
fwest = np.fft.fftfreq(nfft, d=dt/1000)

# Create wavelet in time
wavest0 = np.real(np.fft.ifft(wavest0_fft)[:ntwav])
wavest0 = np.concatenate((np.flipud(wavest0[1:]), wavest0), axis=0)
wavest0 = wavest0 / wavest0.max()

wavest = wavest0.copy()

# Display wavelet
fig, axs = plt.subplots(1, 2, figsize=(20, 5))
fig.suptitle('Statistical wavelet estimate')
axs[0].plot(fwest[:nfft//2], wavest0_fft[:nfft//2], 'k')
axs[0].set_title('Frequency')
axs[1].plot(twav, wav, 'r', label='Wav true')
axs[1].plot(twav, wavest, 'k', label='Wav statistical')
axs[1].legend()
axs[1].set_title('Time');
In [66]:
fd_plus = fd_plus_sparse.copy()
d = d_sparse.copy()
In [67]:
sigmas = np.zeros(nouter)
taus = np.zeros(nouter)
scs = np.zeros((nouter, 2))

d = d.ravel()
f1_alt_inv = np.zeros(2*nt2 * nr)
tau = 0
sigma = np.inf

iouter = 0

while iouter < nouter and sigma > optimalsigma:
    # 1. compute current misfit
    Wavop = Convolve1D(ns*nt2, wavest, wav_c, dims=(nt2, ns), dir=0, dtype='float64')
    Wav1op = Convolve1D(ns*nt2, np.flipud(wavest), wav_c, dims=(nt2, ns), dir=0, dtype='float64')

    Mop = VStack([HStack([Wavop, -1*Wop*Rop]),
                  HStack([-1*Wop*Rollop*R1op, Wav1op])])*BlockDiag([Wop, Wop])

    dw = np.zeros(2*nt2 * nr) if iouter == 0 else Mop * f1_alt_inv.ravel()
    res = d - dw
    sigma = np.linalg.norm(res.ravel())
    sigmas[iouter] = sigma
    
    # 2. update tau
    lamda = - np.linalg.norm(Mop.H * res, ord=np.inf) / sigma
    tau -= (sigma - optimalsigma) / lamda
    taus[iouter] = tau
    
    # 3. invert for focusing functions
    f1_alt_inv = SPGL1(Mop, d, tau=tau, 
                       x0=f1_alt_inv.ravel(), 
                       iter_lim=50 if iouter == nouter-1 else 50, 
                       verbosity=1)[0] 
    f1_alt_inv = f1_alt_inv.reshape(2*nt2, nr)
    dest = Mop * f1_alt_inv.ravel()
    res = d - dest
    res = res.reshape(2*nt2, nr)

    Mf = np.vstack((Mop * np.pad(f1_alt_inv[:nt2].ravel(), (0, nt2*nr), 
                                 mode='constant'),
                    Mop * np.pad(f1_alt_inv[nt2:].ravel(), (nt2*nr, 0), 
                                 mode='constant'))).T

    # 4. scale focusing functions to account for amplitude underestimation
    #sc = solve(Mf.T @ Mf, Mf.T @ d)
    sc = np.linalg.lstsq(Mf, d)[0]
    scs[iouter] = sc
    print(sc)
    
    f1_alt_inv_sc = f1_alt_inv.ravel() * np.hstack((sc[0]*np.ones(nt2*nr), sc[1]*np.ones(nt2*nr)))
    f1_alt_inv_sc = f1_alt_inv_sc.reshape(2*nt2, nr)
    
    dest_sc = Mop * f1_alt_inv_sc.ravel()
    res_sc = d - dest_sc
    res_sc = res_sc.reshape(2*nt2, nr)
    print('res l2', np.linalg.norm(res.ravel()))
    print('res_sc l2', np.linalg.norm(res_sc.ravel()))

    fig, axs = plt.subplots(1, 3, sharey=True, figsize=(15, 9))
    axs[0].imshow(d.reshape(2*nt2, nr), cmap='gray', 
              vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
    axs[0].set_title('data')
    axs[0].axis('tight')
    axs[1].imshow(res.reshape(2*nt2, nr), cmap='gray', 
              vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
    axs[1].set_title('res')
    axs[1].axis('tight')
    axs[2].imshow(res_sc.reshape(2*nt2, nr), cmap='gray', 
              vmin=-clip, vmax=clip, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
    axs[2].set_title('res scaled')
    axs[2].axis('tight')

    f1_alt_inv_tot = f1_alt_inv + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus))
    f1_alt_inv_minus, f1_alt_inv_plus =  f1_alt_inv_tot[:nt2].T, f1_alt_inv_tot[nt2:].T
    
    f1_alt_inv_sc_tot = f1_alt_inv_sc + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus))
    f1_alt_inv_sc_minus, f1_alt_inv_sc_plus =  f1_alt_inv_sc_tot[:nt2].T, f1_alt_inv_sc_tot[nt2:].T
    
    if quickplot:
        f1_alt_inv_tot_plot = np.concatenate((f1_alt_inv_tot[:(2*nt-1)], f1_alt_inv_tot[(2*nt-1):]), axis=1).T
        fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 7))
        ax.imshow(f1_alt_inv_tot_plot.T, cmap='gray', 
                  vmin=-1e5, vmax=1e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
        ax.set_title(r'$f_{inv}$'), axs[0].set_xlabel(r'$x_R$')
        ax.axis('tight')
        ax.set_ylim(1, -1);

    # 5. invert for wavelet
    W1op = Diagonal(w.flatten())
    F1plus = VStack([Convolve1D(nt2, f1_alt_inv_sc[nt2:, ir], offset=nt-2) for ir in range(nr)])
    F1minus = VStack([Convolve1D(nt2, f1_alt_inv_sc[:nt2, ir], offset=nt-1) for ir in range(nr)])
    S = Flip(nt2)
    Fw = VStack([W1op*F1minus, W1op*F1plus*S])

    d1 = Wop*Rop*f1_alt_inv_sc_plus.T.flatten()
    d2 = Wop*R1op*f1_alt_inv_sc_minus.T.flatten()
    dw = np.concatenate((d1.reshape(nt2, ns).T, d2.reshape(nt2, ns).T))

    wavest = np.pad(wavest, (nt-ntwav, nt-ntwav), mode='constant')
    wavest = lsqr(Fw, dw.flatten(), x0=wavest, iter_lim=20, show=False)[0]
    resw = dw.flatten() - Fw @ wavest
    resw = resw.reshape(2*ns, nt2)
    wavest = wavest[nt-ntwav:nt+ntwav-1]
    
    if quickplot:
        plt.figure()
        plt.plot(twav, wav, 'k', lw=2, label='True')
        plt.plot(twav, wavest, 'r', label='Est')
        plt.plot(twav, wavest0, '--b', label='Est it=0')
        plt.legend()
        
        fig, axs = plt.subplots(1, 2, sharey=True, figsize=(5, 7))
        axs[0].imshow(dw.T, cmap='gray', 
                  vmin=-1e5, vmax=1e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
        axs[0].set_title(r'$d_{w}$'), axs[0].set_xlabel(r'$x_R$')
        axs[0].axis('tight')
        axs[0].set_ylim(1, -1);
        axs[1].imshow(resw.T, cmap='gray', 
                  vmin=-1e5, vmax=1e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
        axs[1].set_title(r'$r_{w}$'), axs[0].set_xlabel(r'$x_R$')
        axs[1].axis('tight')
        axs[1].set_ylim(1, -1);

    iouter += 1
================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 1.06e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 1.06e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :       50


ERROR EXIT -- Too many iterations

Products with A     :      19        Total time   (secs) :    45.2
Products with A^H   :      50        Project time (secs) :     4.4
Newton iterations   :       0        Mat-vec time (secs) :    38.1
Line search its     :      18        Subspace iterations :       0
[1.267891 0.      ]
res l2 11359255.772958953
res_sc l2 9901363.193435274

================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 1.76e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 1.76e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :       50


ERROR EXIT -- Too many iterations

Products with A     :      28        Total time   (secs) :    37.2
Products with A^H   :      50        Project time (secs) :     3.8
Newton iterations   :       0        Mat-vec time (secs) :    31.3
Line search its     :      27        Subspace iterations :       0
[1.06583729 1.41790647]
res l2 4225341.608479
res_sc l2 3757611.1693292935

================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 2.26e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 2.26e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :       50


ERROR EXIT -- Too many iterations

Products with A     :      23        Total time   (secs) :    39.2
Products with A^H   :      50        Project time (secs) :     3.9
Newton iterations   :       0        Mat-vec time (secs) :    32.9
Line search its     :      22        Subspace iterations :       0
[1.02220112 1.11915202]
res l2 2158163.38453004
res_sc l2 2038380.022501756

================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 2.55e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 2.55e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :       50


ERROR EXIT -- Too many iterations

Products with A     :      27        Total time   (secs) :    39.6
Products with A^H   :      50        Project time (secs) :     4.0
Newton iterations   :       0        Mat-vec time (secs) :    33.2
Line search its     :      26        Subspace iterations :       0
[1.01056528 1.05324525]
res l2 1370301.6270656495
res_sc l2 1327255.6570137849

================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 2.68e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 2.68e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :       50


ERROR EXIT -- Too many iterations

Products with A     :      30        Total time   (secs) :    37.7
Products with A^H   :      50        Project time (secs) :     3.9
Newton iterations   :       0        Mat-vec time (secs) :    31.6
Line search its     :      29        Subspace iterations :       0
[1.00667948 1.03497246]
res l2 1091638.5942232304
res_sc l2 1068933.4280142237

================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 2.77e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 2.77e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :       50


ERROR EXIT -- Too many iterations

Products with A     :      29        Total time   (secs) :    39.2
Products with A^H   :      50        Project time (secs) :     4.0
Newton iterations   :       0        Mat-vec time (secs) :    33.0
Line search its     :      28        Subspace iterations :       0
[1.00467403 1.0247305 ]
res l2 939020.3589183663
res_sc l2 925911.3645839299

================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 2.83e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 2.83e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :       50


ERROR EXIT -- Too many iterations

Products with A     :      29        Total time   (secs) :    40.4
Products with A^H   :      50        Project time (secs) :     4.2
Newton iterations   :       0        Mat-vec time (secs) :    33.9
Line search its     :      28        Subspace iterations :       0
[1.00346765 1.01821575]
res l2 842121.2448233048
res_sc l2 834139.4492362108

================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 2.89e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 2.89e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :       50


ERROR EXIT -- Too many iterations

Products with A     :      28        Total time   (secs) :    38.4
Products with A^H   :      50        Project time (secs) :     3.9
Newton iterations   :       0        Mat-vec time (secs) :    32.2
Line search its     :      27        Subspace iterations :       0
[1.00268098 1.0140197 ]
res l2 778305.7245568518
res_sc l2 773181.2498350253

================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 2.94e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 2.94e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :       50

Restoring best iterate to objective 731372.468601282

ERROR EXIT -- Too many iterations

Products with A     :      27        Total time   (secs) :    44.4
Products with A^H   :      51        Project time (secs) :     4.6
Newton iterations   :       0        Mat-vec time (secs) :    37.4
Line search its     :      25        Subspace iterations :       0
[1.00212475 1.01189912]
res l2 731372.468601282
res_sc l2 727714.115460439

================================================================================
SPGL1
================================================================================
No. rows              :   322998     
No. columns           :   322998

Initial tau           : 2.98e+08     
Two-norm of b         : 2.81e+07

Optimality tol        : 1.00e-04     
Target one-norm of x  : 2.98e+08

Basis pursuit tol     : 1.00e-06     
Maximum iterations    :       50


ERROR EXIT -- Too many iterations

Products with A     :      26        Total time   (secs) :    38.3
Products with A^H   :      50        Project time (secs) :     4.0
Newton iterations   :       0        Mat-vec time (secs) :    32.0
Line search its     :      25        Subspace iterations :       0
[1.00178489 1.01009692]
res l2 695088.8244207612
res_sc l2 692357.1876396502
In [68]:
g_alt_inv = Gop*f1_alt_inv_tot.flatten()
g_alt_inv = g_alt_inv.reshape(2*nt2, ns)

f1_sparse_inv_minus, f1_sparse_inv_plus =  f1_sparse_inv_tot[:nt2].T, f1_sparse_inv_tot[nt2:].T
g_alt_inv_minus, g_alt_inv_plus =  -g_alt_inv[:nt2].T, np.fliplr(g_alt_inv[nt2:].T)
g_alt_inv_tot = g_alt_inv_minus + g_alt_inv_plus

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_alt_inv_minus.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_alt_inv_plus.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].plot(r[0], directVS,'b'),axs[2].plot(r[0], -directVS,'b')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow((1-w.T)*g_alt_inv_minus.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_alt_inv_plus.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);

fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-clip, vmax=clip, 
           extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_alt_inv_tot.T, cmap='gray', vmin=-clip, vmax=clip, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(t**2*Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(t**2*g_alt_inv_tot[nr//2, nt-1:]/g_alt_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);
In [69]:
fig, axs = plt.subplots(1, 3, figsize=(15, 4))
axs[0].plot(taus, 'k')
axs[0].set_title('tau')
axs[1].plot(scs[:, 0], 'k', label='1st term')
axs[1].plot(scs[:, 1], '--k', label='2nd term')
axs[1].set_title('sc')
axs[2].plot(sigmas, 'k')
axs[2].axhline(optimalsigma, color='k', lw=5)
axs[2].set_title('sigma');

fig, ax = plt.subplots(1, 1, figsize=(10, 3))
ax.plot(taus, sigmas, 'k')
ax.title('tau vs sigma')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-69-3d577d2b0b6c> in <module>
     11 fig, ax = plt.subplots(1, 1, figsize=(10, 3))
     12 ax.plot(taus, sigmas, 'k')
---> 13 ax.title('tau vs sigma')

TypeError: 'Text' object is not callable

iouter=0 tau=0

d = d.ravel() f1_alt_inv = np.zeros(2nt2 nr)

1. compute current misfit

Wavop = Convolve1D(nsnt2, wavest, wav_c, dims=(nt2, ns), dir=0, dtype='float64') Wav1op = Convolve1D(nsnt2, np.flipud(wavest), wav_c, dims=(nt2, ns), dir=0, dtype='float64')

Mop = VStack([HStack([Wavop, -1WopRop]), HStack([-1WopRollopR1op, Wav1op])])BlockDiag([Wop, Wop])

dw = np.zeros(2nt2 nr) if iouter == 0 else Mop * f1_alt_inv.ravel() res = d - dw sigma = np.linalg.norm(res.ravel())

2. update tau

lamda = - np.linalg.norm(Mop.H * res, ord=np.inf) / sigma tau -= (sigma - optimalsigma) / lamda

3. invert for focusing functions

f1_alt_inv = SPGL1(Mop, d, tau=tau, iter_lim=10, verbosity=1)[0] # x0=f1_alt_inv.ravel(), f1_alt_inv = f1_alt_inv.reshape(2nt2, nr) dest = Mop f1_alt_inv.ravel()

Mf = np.vstack((Mop np.pad(f1_alt_inv[:nt2].ravel(), (0, nt2nr), mode='constant'), Mop np.pad(f1_alt_inv[nt2:].ravel(), (nt2nr, 0), mode='constant'))).T

4. scale focusing functions to account for amplitude underestimation

sc = solve(Mf.T @ Mf, Mf.T @ d)

sc = np.linalg.lstsq(Mf, d)[0] print(sc) f1_alt_inv = f1_alt_inv.ravel() np.hstack((sc[0]np.ones(nt2nr), sc[1]np.ones(nt2nr))) f1_alt_inv = f1_alt_inv.reshape(2nt2, nr)

dest_sc = Mop * f1_alt_inv.ravel()

ffig, axs = plt.subplots(1, 3, sharey=True, figsize=(15, 9)) axs[0].imshow(d.reshape(2nt2, nr), cmap='gray', vmin=-1e6, vmax=1e6, extent=(r[0,0], 2r[0,-1], t[-1], -t[-1])) axs[0].axis('tight') axs[1].imshow(dest.reshape(2nt2, nr), cmap='gray', vmin=-1e6, vmax=1e6, extent=(r[0,0], 2r[0,-1], t[-1], -t[-1])) axs[1].axis('tight') axs[2].imshow(dest_sc.reshape(2nt2, nr), cmap='gray', vmin=-1e6, vmax=1e6, extent=(r[0,0], 2r[0,-1], t[-1], -t[-1])) axs[2].axis('tight')

f1_alt_inv_tot = f1_alt_inv + np.concatenate((np.zeros((2*nt-1, nr)), fd_plus)) f1_alt_inv_minus, f1_alt_inv_plus = f1_alt_inv_tot[:nt2].T, f1_alt_inv_tot[nt2:].T

f1_alt_inv_tot_plot = np.concatenate((f1_alt_inv_tot[:(2nt-1)], f1_alt_inv_tot[(2nt-1):]), axis=1).T fig, ax = plt.subplots(1, 1, sharey=True, figsize=(5, 7)) ax.imshow(f1_alt_inv_tot_plot.T, cmap='gray', vmin=-1e5, vmax=1e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1])) ax.settitle(r'$f{inv}$'), axs[0].set_xlabel(r'$x_R$') ax.axis('tight') ax.set_ylim(1, -1);

5. invert for wavelet

W1op = Diagonal(w.flatten()) F1plus = VStack([Convolve1D(nt2, f1_alt_inv[nt2:, ir], offset=nt-2) for ir in range(nr)]) F1minus = VStack([Convolve1D(nt2, f1_alt_inv[:nt2, ir], offset=nt-1) for ir in range(nr)]) S = Flip(nt2) Fw = VStack([W1opF1minus, W1opF1plus*S])

d1 = WopRopf1_alt_inv_plus.T.flatten() d2 = WopR1opf1_alt_inv_minus.T.flatten() dw = np.concatenate((d1.reshape(nt2, ns).T, d2.reshape(nt2, ns).T))

wavest = np.pad(wavest, (nt-ntwav, nt-ntwav), mode='constant') wavest = lsqr(Fw, dw.flatten(), iter_lim=5, show=True)[0] # x0=wavest wavest = wavest[nt-ntwav:nt+ntwav-1]

if quickplot: plt.figure() plt.plot(twav, wav, 'k', lw=2) plt.plot(twav, wavest, '--r');

6. remove scaling from focusing functions

f1_alt_inv = f1_alt_inv.ravel() np.hstack((sc[0]np.ones(nt2nr), sc[1]np.ones(nt2nr))) f1_alt_inv = f1_alt_inv.reshape(2nt2, nr)

In [ ]: